%This file aims to create IVs---averaging market share and market dominance
%indicators across all outside associated markets accounting for hospital 
%patient size.


clear
tic

textFileName1=sprintf('sysyr.raw'); %the input file: the hsayr_id and sysyr_id and numhsa for involved multi-region system
textFileName2=sprintf('input4createIV.raw'); %the input file of basic info of all hospitals
sysyr=dlmread(textFileName1);  %the matrix: hsayr_id sysyr_id numhsa
H=dlmread(textFileName2);  %the matrix: sysyr_id ahayr_id hsayr_id prechs adj_vnd bdtot sysid nprofit profit teaching ncompet ratio perc_mcr perc_mcd
choice=1:13;
Tchoice=choice';                     
ns=numel(choice);



%%%%%%%%%%%%%%%%%%%%
%The first way is to find out the largest HSA according to how many
%member hospitals locate in.
exo_mktdom=zeros(size(sysyr,1),ns);
exo_mktsh=zeros(size(sysyr,1),ns);
for i=1:size(sysyr,1)
    tempsys=H(H(:,1)==sysyr(i,2),:);
    %%%%%%%%%%
    %The following is to compute the exogenous big hsa for a hospital
    %system
    syshsa=sort(tempsys(:,3));
    unique_syshsa=unique(syshsa);
    occur=histc(syshsa,unique_syshsa);
    bighsa=unique_syshsa(occur==max(occur));
    exo_bighsa=bighsa;
    %exo_bighsa(exo_bighsa(:)==sysyr(i,1))=[];
    if min(abs(exo_bighsa-sysyr(i,1)))==0
        exo_bighsa=[];
    end
  if isempty(exo_bighsa)==0 
      part_mktdom=zeros(numel(exo_bighsa),ns);
      part_mktsh=zeros(numel(exo_bighsa),ns);
      for m=1:numel(exo_bighsa)
          prechs=H(H(:,3)==exo_bighsa(m),4);
          prechs(prechs==1)=[];
          if isempty(prechs)==0
          occur=histc(prechs,choice);
          occur=reshape(occur,1,[]);           
          part_mktsh(m,:)=occur/sum(occur); %a vector 1 by 13
          temp_prechs=prechs;
          temp_prechs(prechs==2|prechs==13)=[];
          if isempty(temp_prechs)==0
             occur=histc(temp_prechs,choice)';
             mktdom=choice(occur==max(occur));
             vmktff=max((~(Tchoice(:,ones(size(mktdom,2),1))-mktdom(ones(size(choice,2),1),:))),[],2);
             part_mktdom(m,:)=vmktff';            
          end 
          end
      end
      exo_mktdom(i,:)=mean(part_mktdom,1);
      exo_mktsh(i,:)=mean(part_mktsh,1);
  end 
end


iv_mktdom=[sysyr(:,1:2),exo_mktdom];
iv_mktdom(~any(iv_mktdom(:,3:end),2),:)=[];  %### UPDATED on 07/25/2019
iv_mktsh=[sysyr(:,1:2),exo_mktsh];
iv_mktsh(~any(iv_mktsh(:,3:end),2),:)=[];  %### UPDATED on 07/25/2019

iv_mktdom=sortrows(iv_mktdom,[1,2]);
iv_mktsh=sortrows(iv_mktsh,[1,2]);

hsayr=unique(sort(iv_mktdom(:,1)));

output_ivdom=zeros(numel(hsayr),ns);
output_ivmsh=zeros(numel(hsayr),ns);
for i=1:numel(hsayr)
    output_ivdom(i,:)=mean(iv_mktdom(iv_mktdom(:,1)==hsayr(i),3:end),1);
    output_ivmsh(i,:)=mean(iv_mktsh(iv_mktsh(:,1)==hsayr(i),3:end),1);   
end


textFileName11=sprintf('iv_dom_tmp.csv');
dlmwrite(textFileName11,[output_ivdom,hsayr],'delimiter', ',', 'precision', 8)
textFileName12=sprintf('iv_msh_tmp.csv');
dlmwrite(textFileName12,[output_ivmsh,hsayr],'delimiter', ',', 'precision', 8)
%}




%%%%%%%%%%%%%%%%%%%% ### UPDATED on 11/16/2020: construct biggest HSA and
%%%%%%%%%%%%%%%%%%%% IV, BOTH based on bedsize-adjusted mktsh
exo_mktdombed=zeros(size(sysyr,1),ns);
exo_mktshbed=zeros(size(sysyr,1),ns);
for i=1:size(sysyr,1)
    tempsys=H(H(:,1)==sysyr(i,2),:);
    tempsys=sortrows(tempsys,[3,6]);
    %The following is to compute the exogenous big hsa for a hospital
    %system
    hsabed=tempsys(:,[3,6]);  %matrix: hsayr_id bdtot
    hsabed=sortrows(hsabed,[1,2]);
    [hsa,~,idhsa]=unique(hsabed(:,1));  %matrix: unique_hsa id_in_unique_hsa_for_hsabed
    accumbed=accumarray(idhsa,hsabed(:,2));
    bighsa=hsa(accumbed==max(accumbed));
    exo_bighsa=bighsa;
    %exo_bighsa(exo_bighsa(:)==sysyr(i,1))=[];
    if min(abs(exo_bighsa-sysyr(i,1)))==0
        exo_bighsa=[];
    end
  if isempty(exo_bighsa)==0 
      part_mktdom=zeros(numel(exo_bighsa),ns);
      part_mktsh=zeros(numel(exo_bighsa),ns);
      for m=1:numel(exo_bighsa)
          fulmsh=[choice',zeros(ns,1)];
          tempprechs=H(H(:,3)==exo_bighsa(m),[4,6]);  %matrix: prechs bdtot
          tempprechs=sortrows(tempprechs,[1,2]);
          tempprechs(tempprechs(:,1)==1,:)=[];
          if isempty(tempprechs)==0
              [unqprechs,~,idprechs]=unique(tempprechs(:,1));
              accumprchsbd=accumarray(idprechs,tempprechs(:,2));
              tmpmsh=accumprchsbd/sum(accumprchsbd);
              idx=sub2ind(size(fulmsh),unqprechs,2*ones(numel(unqprechs),1));
              fulmsh(idx)=tmpmsh;
              part_mktsh(m,:)=fulmsh(:,2);
              fulmsh([2,13],2)=0;
              if sum(fulmsh(:,2))>0               
                 mktdom=fulmsh(fulmsh(:,2)==max(fulmsh(:,2)),1);
                 vmktff=max((~(Tchoice(:,ones(size(mktdom,2),1))-mktdom(ones(size(choice,2),1),:))),[],2);
                 part_mktdom(m,:)=vmktff';            
              end 
          end
      end
      exo_mktdombed(i,:)=mean(part_mktdom,1);
      exo_mktshbed(i,:)=mean(part_mktsh,1);
  end 
end
    
iv_mktdombed=[sysyr(:,1:2),exo_mktdombed];
iv_mktdombed(~any(iv_mktdombed(:,3:end),2),:)=[];  %### UPDATED on 07/25/2019
iv_mktshbed=[sysyr(:,1:2),exo_mktshbed];
iv_mktshbed(~any(iv_mktshbed(:,3:end),2),:)=[];  %### UPDATED on 07/25/2019

iv_mktdombed=sortrows(iv_mktdombed,[1,2]);
iv_mktshbed=sortrows(iv_mktshbed,[1,2]);

hsayr=unique(sort(iv_mktdombed(:,1)));

output_ivdombed=zeros(numel(hsayr),ns);
output_ivmshbed=zeros(numel(hsayr),ns);
for i=1:numel(hsayr)
    output_ivdombed(i,:)=mean(iv_mktdombed(iv_mktdombed(:,1)==hsayr(i),3:end),1);
    output_ivmshbed(i,:)=mean(iv_mktshbed(iv_mktshbed(:,1)==hsayr(i),3:end),1);   
end


textFileName11bedtmp=sprintf('iv_dombed_tmp.csv');
dlmwrite(textFileName11bedtmp,[output_ivdombed,hsayr],'delimiter', ',', 'precision', 8)
textFileName12bedtmp=sprintf('iv_mshbed_tmp.csv');
dlmwrite(textFileName12bedtmp,[output_ivmshbed,hsayr],'delimiter', ',', 'precision', 8)
%}





%%%%%%%%%%%%%%%%%%%% ### UPDATED on 11/16/2020: construct biggest HSA and
%%%%%%%%%%%%%%%%%%%% IV, ONLY IV based on bedsize-adjusted mktsh
exo_mktdomIVbed=zeros(size(sysyr,1),ns);
exo_mktshIVbed=zeros(size(sysyr,1),ns);
for i=1:size(sysyr,1)
    tempsys=H(H(:,1)==sysyr(i,2),:);
    %%%%%%%%%%
    %The following is to compute the exogenous big hsa for a hospital
    %system
    syshsa=sort(tempsys(:,3));
    unique_syshsa=unique(syshsa);
    occur=histc(syshsa,unique_syshsa);
    bighsa=unique_syshsa(occur==max(occur));
    exo_bighsa=bighsa;
    %exo_bighsa(exo_bighsa(:)==sysyr(i,1))=[];
    if min(abs(exo_bighsa-sysyr(i,1)))==0
        exo_bighsa=[];
    end
  if isempty(exo_bighsa)==0 
      part_mktdom=zeros(numel(exo_bighsa),ns);
      part_mktsh=zeros(numel(exo_bighsa),ns);
      for m=1:numel(exo_bighsa)
          fulmsh=[choice',zeros(ns,1)];
          tempprechs=H(H(:,3)==exo_bighsa(m),[4,6]);  %matrix: prechs bdtot
          tempprechs=sortrows(tempprechs,[1,2]);
          tempprechs(tempprechs(:,1)==1,:)=[];
          if isempty(tempprechs)==0
              [unqprechs,~,idprechs]=unique(tempprechs(:,1));
              accumprchsbd=accumarray(idprechs,tempprechs(:,2));
              tmpmsh=accumprchsbd/sum(accumprchsbd);
              idx=sub2ind(size(fulmsh),unqprechs,2*ones(numel(unqprechs),1));
              fulmsh(idx)=tmpmsh;
              part_mktsh(m,:)=fulmsh(:,2);
              fulmsh([2,13],2)=0;
              if sum(fulmsh(:,2))>0               
                 mktdom=fulmsh(fulmsh(:,2)==max(fulmsh(:,2)),1);
                 vmktff=max((~(Tchoice(:,ones(size(mktdom,2),1))-mktdom(ones(size(choice,2),1),:))),[],2);
                 part_mktdom(m,:)=vmktff';            
              end 
          end
      end
      exo_mktdomIVbed(i,:)=mean(part_mktdom,1);
      exo_mktshIVbed(i,:)=mean(part_mktsh,1);
  end 
end
    
iv_mktdomIVbed=[sysyr(:,1:2),exo_mktdomIVbed];
iv_mktdomIVbed(~any(iv_mktdomIVbed(:,3:end),2),:)=[];  %### UPDATED on 07/25/2019
iv_mktshIVbed=[sysyr(:,1:2),exo_mktshIVbed];
iv_mktshIVbed(~any(iv_mktshIVbed(:,3:end),2),:)=[];  %### UPDATED on 07/25/2019

iv_mktdomIVbed=sortrows(iv_mktdomIVbed,[1,2]);
iv_mktshIVbed=sortrows(iv_mktshIVbed,[1,2]);

hsayr=unique(sort(iv_mktdomIVbed(:,1)));

output_ivdomIVbed=zeros(numel(hsayr),ns);
output_ivmshIVbed=zeros(numel(hsayr),ns);
for i=1:numel(hsayr)
    output_ivdomIVbed(i,:)=mean(iv_mktdomIVbed(iv_mktdomIVbed(:,1)==hsayr(i),3:end),1);
    output_ivmshIVbed(i,:)=mean(iv_mktshIVbed(iv_mktshIVbed(:,1)==hsayr(i),3:end),1);   
end


textFileName11IVbed=sprintf('iv_dombed.csv');
dlmwrite(textFileName11IVbed,[output_ivdomIVbed,hsayr],'delimiter', ',', 'precision', 8)
textFileName12IVbed=sprintf('iv_mshbed.csv');
dlmwrite(textFileName12IVbed,[output_ivmshIVbed,hsayr],'delimiter', ',', 'precision', 8)
%}





toc









